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ABSTRACT 


Research  under  this  contract  has  been  concentrated  on  major  problems 
in  numerical  linear  algebra,  (i)  The  determination  of  error  bounds  for 
Gaussian  elimination,  (ii)  The  generalized  eigenvalue  problem  Ax  *  ABx 
and  its  natural  extension  to  the  computation  of  the  canonical  form  of  the 
pencil  A  -  XB  where  A  and  B  are  mxn  matrices.  In  addition,  the  numerical 
aspects  of  various  problems  in  linear  system  theory  and  related  fields  have 
been  studied. 

RESEARCH  REPORT 

1.  A  Posteriori  Error  Bounds  for  Gaussian  Elimination 

Although  in  principle  it  has  been  known  for  many  years  how  to  derive  error 
bounds  for  the  solution  of  a  system  of  linear  equations  calculated  by  Gaussian 
elimination,  programs  in  this  area  have  mainly  been  in  machine  code  and  on  an 
ad  hoc  basis.  Collaboration  with  F.  W.  J.  Olver  of  the  Uiversity  of  Maryland 
has  led  to  the  development  of  programs  which  can  be  coded  in  high  level 
languages  and  should  prove  fully  portable.  They  require  no  special  hardware 
or  software.  The  steps  may  be  described  in  matrix  terms  by  the  following 
equations. 

MA^  *  A  ,  Mb^  *  b  .  (Gaussian  elimination) 

M  unit  lower  triangular;  A^  upper  triangular 

VA(2n-i)  _  A(n)  ^  vb(2n-l)  _  fe(n)  . 

Here  V  is  unit  upper  triangular  and  A^n  ^  is  diagonal. 
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i.e.  Che  solution  of  a  diagonal  system. 

A  backward  error  analysis  is  performed  on  each  of  these  steps  giving 
relations  of  the  type 

MA(n)  -  A  -  AA  ,  Mb(n)  '  b  -  Ab 

where  M  ,  A^  are  the  computed  matrices  and  AA  and  Ab  cover  the  rounding 
errors.  For  |AA|  etc.  we  have  only  error  bounds.  For  the  determination  of 
the  error  bounds  we  require  the  solution  of  the  equations 
LM  =  I  and  UV  «  I 

and  the  error  analysis  of  these  processes  gives  for  the  computed  matrices 
L  and  U 

LM-I+AL  ,  UV-I+AU  . 

The  use  of  these  equations  enables  a  bound  to  be  found  for  x  -  x  . 

The  inherent  errors  in  the  data  may  be  included  in  AA  and  Ab.  The  analysis 
provides  for  the  possible  use  of  lower  precision  arithmetic  in  the  determination 
of  L  and  U  and  in  the  computation  of  the  error  bounds  themselves.  This  work 
is  the  subject  of  a  joint  paper  (1)  to  be  submitted  for  publication. 

The  collaboration  with  Olver,  who  had  not  previously  worked  in  the 
linear  algebra  field,  has  been  very  productive  for  me.  It  came  at  a  time  when 
1  was  already  interested  in  the  problem  of  making  research  on  rounding  error 
analysis  more  readily  accessible  to  non-specialists  in  the  field,  partly 
as  a  result  of  giving  lectures  in  the  Statistical-Numerical  Analysis  Summer 
Schools  at  Delaware  (1980  and  1981).  I  am  now  convinced  that  rounding  error 
analysis  should  be  presented  in  an  entirely  different  way.  It  was  impractical 
to  let  this  new  approach  influence  the  presentation  in  the  joint  paper  with 
Olver  but  it  will  be  the  main  thesis  of  a  monograph  on  Rounding  Error  Analysis 
which  is  now  in  preparation.  It  is  interesting  that  in  many  areas  the  whole 


theory  can  be  illustrated  by  means  of  examples  with  small  Integers  in  which 
the  only  errors  are  integer  errors,  so  that  no  true  rounding  error  analysis 
is  involved.  Hitherto  the  rather  tedious  detail  associated  with  rounding 
errors  has  tended  to  obscure  the  underlying  simplicity  of  what  is  being 
done  and  the  underlying  mathematical  analysis  is  submerged. 

J.  H.  Wilkinson 

2.  The  Calculation  of  Error  Bounds  for  Computed  Eigenvalues  and  Eigenvectors 
and  Invariant  Subspaces 

In  practice  it  is  rare  for  the  complete  eigensystem  of  a  large  matrix 
to  be  required.  Commonly  attention  is  focussed  on  a  few  key  eigenvalues  and 
eigenvectors.  An  algorithm  is  therefore  desirable  for  deriving  error  bounds 
for  a  single  eigenpair  (i.e.  value  and  vector)  without  requiring  information 
on  the  remainder  of  the  system.  A  method  for  deriving  rigorous  error  bounds 
for  such  a  pair  X  and  x  was  developed  and  has  been  published  in  Numerische 
Mathematik  (2) .  A  pleasing  feature  is  that  in  the  process  of  deriving  the  bounds 
an  improved  eigenpair  is  determined  and  the  bound  is  for  this  eigenpair.  The 
method  covers  both  Ax  »  Xx  and  Ax  *  XBx  .  In  addition  to  a  completely  rigorous 
method  a  more  practical  technique  is  described  which  gives  realistic  bounds 
with  far  less  computation. 

A  less  pleasing  feature  is  that  the  bound  (and  indeed,  the  true  error)  is 
limited  by  the  condition  number  of  a  certain  matrix  C.  The  matrix  C  is  ill- 
conditioned  when  X  is  close  to  multiple  roots.  The  latter  are  not  necessarily 
Ill-conditioned  and  hence  this  is  an  inherent  weakness  of  the  method. 
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This  weakness  has  been  overcome  by  determining  generators  of  the  invariant 
subspace  associated  with  clusters  of  eigenvalues.  This  material  was  presented 
as  an  invited  paper  at  an  international  symposium  in  honor  of  H.  Rutishauser 
at  ETH  Zurich  and  has  been  published  in  the  proceedings  of  the  Symposium  (3). 

Although  the  algorithms  described  in  the  above  two  papers  were  essentially 
iterative  they  were  designed  to  deal  with  initial  approximations  of  high 
accuracy  such  as  are  provided  by  stable  algorithms  for  solving  the  eigenvalue 
problem.  Hence  it  was  envisaged  that  only  one  of  two  iterations  were  required 
to  obtain  the  limiting  accuracy.  Although  they  could  be  used  for  complex 
eigenvalues  simply  by  working  in  complex  arithmetic  they  were  inefficient  in 
dealing  with  complex  conjugate  pairs  of  eigenvalues  of  real  matrices.  In  a 
third  paper  (4)  resulting  from  collaboration  with  J.  Dongarra  of  Argonne 
National  Laboratory  and  C.  Moler  of  the  University  of  New  Mexico  extensions 
of  the  basic  algorighm  are  described.  These  give  (i)  Quadratic  convergence, 
efficiency  being  achieved  in  the  iterations  by  an  updating  technique,  (ii) 

Full  efficienty  for  complex  conjugate  eigenvalues,  (iii)  Accurate  generators 
of  invariant  subspaces  associated  with  close  clusters  of  eigenvalues  whether  or 
they  correspond  to  linear  elementary  divisors.  Analyses  of  the  influence  of 
rounding  errors  on  the  performance  of  quadratically  convergent  algorithms  are 
often  based  on  unrealistic  assumptions,  mainly  because  of  the  tedious  nature  of 
the  algebra  involved.  In  a  final  section  of  this  paper  a  program  is  described 
for  performing  this  analysis  on  a  computer,  thereby  making  it  a  practical 
proposition  to  use  a  rigorous  rounding  error  analysis. 


J.  H.  Wilkinson 


3.  The  Determination  of  the  Distance  of  a  Matrix  from  the  Nearest  Defective 
Matrix  and  the  Corresponding  Problem  for  Ax  -  XBx. 

This  is  a  problem  of  considerable  interest  to  control  engineers.  Early 
work  in  connection  with  the  standard  problem  was  done  by  Kahan,  Ruhe,  and 
Wilkinson.  The  techniques  developed  by  them  have  been  greatly  Improved  and 
extended  to  the  generalized  problem  (the  more  important  case).  In  the  course 
of  this  work  rather  general  results  were  derived  which  should  be  of  considerable 
value  in  the  backward  error  analysis  of  eigenvalue  algorithms.  This  work  was 
presented  as  an  invited  lecture  at  an  International  meeting  held  at  the  Univerity 
of  Manitoba  and  the  paper  (5)  has  been  published  in  the  proceedings. 

Work  has  continued  on  problems  associated  with  Knonecker's  canonical 
from  arising  in  control  theory.  Backward  stable  algor ighms  have  been  derived 
for  computing  Knonecker's  canonical  form  and  also  the  Drazin  Inverse.  The 
techniques  used  in  these  papers  are  now  being  widely  adopted.  A  summary  of 
this  work  was  presented  at  the  8th  Gatlinburg  meeting  on  Numerical  Linear 
Algebra  held  at  Oxford  in  July  1981;  a  full  session  at  this  meeting  was  devoted 
to  the  numerical  problems  arising  in  the  control  field. 

J.  H.  Wilkinson 

A.  Application  of  the  generalized  Eigenstructure  Problem  in  Linear  System 

Theory 

A  new  algorithmic  approach  was  developed  recently  to  tackle  the  generalized 
eigenstructure  problem  in  its  most  general  form.  We  applied  these  ideas  to  a  set 
of  specific  problems.  For  each  of  them  a  numerical  analysis  of  the  problem 
was  performed  (6),  (7),  (8),  (9),  (10).  Some  algorithms  were  also  implemented 
and  tested  (9),  (10).  The  specific  problems  treated  were: 
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-  Computation  of  the  controllable  and  unobservable  subspaces,  and 
of  the  Kalman  decomposition. 

-  Computation  of  the  supremal  (A, B) -invariant  subspace  and  (A,B)- 
controllability  subspace  in  the  kernel  of  C,  computation  of  zeros 
(7), (9), (12). 

-  Construction  and  analysis  of  generalized  state-space  systems  (7). 

-  Computing  zeros  of  polynominal  system  models  (6) . 

-  Construction  of  cascade  and  spectral  factorizations  (8), (10). 

-  Solving  Riccati  equations  arising  in  linear  system  theory  (10). 

Paul  Van  Dooren 

5.  Structural  Properties  and  Sensitivity  of  Balanced  Realizations 

The  robustness  of  certain  properties  (such  as  poles,  zeros,  observability, 
controllablity,  step  response,  etc.)  of  systems  in  state-space  form  depends 
highly  on  the  chosen  coordinate  system.  We  analyzed  the  sensitivity  of  these 
properties  for  the  so-called  "balanced  realization".  The  structural  properties 
of  balanced  realizations  for  J-unitary  transfer  functions  was  investigated  (8) . 
Extensions  for  time-varying  systems  were  considered  (11).  Possible  implications 
for  the  synthesis  of  filters  and  predictors  for  "almost  stationary"  stochastic 
processes  are  under  current  investigations. 


Paul  Van  Dooren 
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